function D = Ising(T)    % T is the tempure.  +-+-+   Ising area.
    nTrials = 100000;    % for stat.          |3|4|   N = 2.
    startup = 100000;    % for steady.        +-+-+   free boundary.
    d = zeros(1, nTrials);   % results.       |1|2| 
    beta = 1/T;              %                +-+-+ 
    M = 0;
    s = 2 * round(rand(1, 4)) - 1;   % randomly set 1 or -1.
    Es = s(1) * s(2) + s(3) * s(4);        
    Es = Es + s(1) * s(3) + s(2) * s(4);   
    Es = -Es;                              
    for t = 1 : (startup + nTrials)        
        k = fix(1 + 4 * rand);  % randomly pick 1 ~ 4.
        y = s;         % y is the suggestion dist in s' neighbour. 
        y(k) = -s(k);  % randomly flip the status on one point.
        Ey = y(1) * y(2) + y(3) * y(4);
        Ey = Ey + y(1) * y(3) + y(2) * y(4);
        Ey = -Ey;
        h = min(1, exp(-beta * (Ey - Es)))
        if (rand < h)
            s = y;
            Es = Ey;
        end
        if (t > startup)
            Ms = s(1) + s(2) + s(3) + s(4);  
            M = M + Ms;                      
            d(t - startup) = Ms;
        end
    end
    x = -4:1:4;
    hist(d, x);
    M = (M / nTrials) / 4;
    D = sum(d.^2);
end